C     ==================================================================       
C     ------Reynolds number
      REAL*8 FUNCTION DRE(G, L, MU)
      REAL*8 G, L, MU
C      
      DRE = G*L/MU
C      
      RETURN
      END
      

C     ==================================================================       
C     ------Heat transfer Nusselt number
      REAL*8 FUNCTION DNU(GR, RE, PR, WVCF, OPT)
      REAL*8 GR, RE, PR, WVCF
      INTEGER OPT
C      
C     If OPT = 1, use tube-side correlation.
C     Otherwise, use shell-side correlation.
      IF (OPT .EQ. 1) THEN
!        IF ((GR .GT. 1.D-1) .AND. (GR .LE. 1.D2)) THEN
!          DNU = 3.66+(0.19*GR**0.8)/(1+0.117*GR**0.467)*WVCF
!        ELSE IF (GR .GT. 1.D2) THEN
!          DNU = 1.86*GR**(1.0/3.0)*WVCF
!        ELSE
          DNU = 3.66+(0.19*GR**0.8)/(1+0.117*GR**0.467)*WVCF
C         WRITE(*, *) 'Tube-side Reynolds number too low'      
!        END IF
      ELSE
!        IF ((RE .GT. 2.D2) .AND. (RE .LT. 2.D5)) THEN
!          DNU = 0.16*RE**0.6*PR**0.33*WVCF**0.14
!        ELSE IF ((RE .GT. 1.0) .AND. (RE .LE. 1.D5)) THEN
!          DNU = (0.4*RE**0.5+0.06*RE**(2.0/3.0))*PR**0.4*WVCF**0.25
!        ELSE
          DNU = (0.4*RE**0.5+0.06*RE**(2.0/3.0))*PR**0.4*WVCF**0.25
C         WRITE(*, *) 'Shell-side Reynolds number too low'      
!        END IF          
      END IF
C
      RETURN
      END


C     ==================================================================       
C     ------Prandtl number
      REAL*8 FUNCTION DPR(CP, MU, K)
      REAL*8 CP, MU, K
C     
      DPR = CP*MU/K
C      
      RETURN
      END
      

C     ==================================================================       
C     ------Graetz number
      REAL*8 FUNCTION DGR(RE, PR, D, L)
      REAL*8 RE, PR, D, L
C      
      DGR = RE*PR*(D/L)
C      
      RETURN
      END
      

C     ==================================================================       
C     ------Latent heat
      REAL*8 FUNCTION DHVAP(T)
      REAL*8 T, TR
      REAL*8 C(4), TC, MW
C     Latent heat of water [J/kg]
C     Input parameter:
C     T - temperature [K] if T < 273.15, T = T+273.15
C      
      C(1) = 5.2053E7
      C(2) = 0.3199
      C(3) = -0.212
      C(4) = 0.25795
      TC = 647.13
      MW = 18.015
C
      TR = T/TC
C
      IF (T .LE. 273.15) THEN
        T = 273.15+T
      END IF
      DHVAP = C(1)*(1-TR)**(C(2)+C(3)*TR+C4*TR*TR)/MW
C
      RETURN
      END

      
C     ==================================================================      
C     ------Gradient of vapor pressure
      REAL*8 FUNCTION GRADSP(T)
      REAL*8 T, DEXP, DLOG
      REAL*8 A, B, C
      REAL*8 C1, C2, C3, C4, C5
C     Gradient of vapor pressure [Pa]
C     Input parameter:
C     T - temperature [K] if T < 273.15, T = T+273.15
C           
      A = 23.238
      B = 3841.0
      C = -45.0
      
      C1 = 73.649
      C2 = -7258.2
      C3 = -7.3037
      C4 = 4.1653E-6
      C5 = 2.0
C
      IF (T .LE. 273.15) THEN
        T = 273.15+T
      END IF
C     Antonie correlation
C     GRADSP = (B/T**2.0)*DEXP(A-B/(C+T))
C
      GRADSP = DEXP(C1+C2/T+C3*DLOG(T)+C4*T**C5)*
     +         (-C2/T**2+C3/T+C4*C5*T**(C5-1.0))      
      RETURN
      END


C     ==================================================================      
      REAL*8 FUNCTION SVP(T)
      REAL*8 T, DEXP, DLOG
      REAL*8 A, B, C
      REAL*8 C1, C2, C3, C4, C5
C     Saturated vapor pressure [Pa]
C     Input parameter:
C     T - temperature [K] if T < 273.15, T = T+273.15
C           
      A = 23.238
      B = 3841.0
      C = -45.0
      
      C1 = 73.649
      C2 = -7258.2
      C3 = -7.3037
      C4 = 4.1653E-6
      C5 = 2.0
C
      IF (T .LE. 273.15) THEN
        T = 273.15+T
      END IF
C     Antonie correlation
C     SVP = DEXP(A-B/(C+T))
C
      SVP = DEXP(C1+C2/T+C3*DLOG(T)+C4*T**C5)
            
      RETURN
      
      END
      
C     ==================================================================      
C     ------Calculate pressure drop
      REAL*8 FUNCTION PDROP(G, L, D, RHO, MU, WVCF, OPT)
C     Input parameters
C     G - flowing mass flux [kg/m2-s]
C     L - total length [m]
C     D - characteristic diameter [m]
C     RHO - fluid density [kg/m3]
C     MU - fluid viscosity [Pa-s]
C     WVCF  - wall viscosity correction factor
C     OPT - tube-side or shell-side switch (1 - tube side)
      REAL*8 G, L, D, RHO, MU, WVCF
      INTEGER OPT
      REAL*8 RE, F, DPT, DPR
C      
      RE = G*D/MU
      
      IF (RE .LT. 2000) THEN
        F = 16.0/RE
      ELSE
        F = 0.0035+0.264/RE**0.41
      END IF
      DPT = 2.0*F*G**2.0/RHO*(L/D)*WVCF**(-0.14)
      IF (OPT .EQ. 1) THEN
        DPR = 2.0*G**2.0/RHO
      ELSE
        DPR = 0.0
      END IF
      PDROP = DPT+DPR
C
      RETURN
      END
      
      
C     ==================================================================  
      SUBROUTINE NumStream(SIDENUM, IONUM, NUMS)
C     Input parameters
C     SIDENUM - 1 for tube side; 2 for shell side
C     IONUM - 1 for inlet; 2 for outlet
      INTEGER SIDENUM, IONUM
C     Output parameter
      INTEGER NUMS
      
      IF (SIDENUM .EQ. 1) THEN ! Tube side
        IF (IONUM .EQ. 1) THEN ! Inlet
          NUMS = 1
        ELSE                   ! Outlet
          NUMS = 2
        END IF
      ELSE                     ! Shell side
        IF (IONUM .EQ. 1) THEN ! Inlet
          NUMS = 2
        ELSE                   ! Outlet
          NUMS = 1
        END IF
      END IF
      
      END
      
      
C     ==================================================================       
      REAL*8 FUNCTION GradDP(d, rho, v, mu)
C     Calculate the gradient of pressure      
      REAL*8 d, rho, v, mu
      REAL*8 Re, f
      
      Re = d*v*rho/mu
      
      IF (Re .LT. 2.1d3) THEN
        f = 16.d0/Re
      ELSE IF (Re .LT. 1.0d5) THEN
        f = 0.0791d0/Re**0.25d0
      ELSE
        WRITE(*, *) 'Reynolds number is out of range in GradDP'
        STOP
      END IF
      
      GradDP = -4.d0/d*(0.5d0*rho*v**2.0d0)*f
      
      END

C     ==================================================================       
      REAL*8 FUNCTION SVMU(T)
C     Correlation obtained from a curve fit of data from 'Heat Transfer'
C     by Alan J. Chapman, 1974
C     T - (input) given temperature [K or C]
C     SVMU - (output) viscosity [kg/m-s]
      REAL*8 T
      REAL*8 P_KPA, P_PSI
      REAL*8 C1, C2, C3, C4
      REAL*8 SVP
      
      DATA C1 /3.14D-2/, C2 /2.9675D-5/, C3 /-1.60583D-8/, 
     &     C4 /3.768986D-12/
      
      P_KPA = SVP(T)/1.D3
      P_PSI = P_KPA/6.894757D0
      SVMU = C1+C2*P_PSI+C3*P_PSI**2.D0+C4*P_PSI**3.D0
      SVMU = SVMU*4.1338D-4
      
      END FUNCTION


C     ==================================================================       
      REAL*8 FUNCTION SVLAMBDA(T)
C     Correlate the thermal conductivity of saturated stream
C     Curve fit from data in "Heat Transfer" by Alan J. Chapman, 1974
C     T - (input) given temperature [K or C]
C     SVLAMBDA - (output) thermal conductivity [W/m-K]
      REAL*8 T
      REAL*8 C1, C2, C3
      REAL*8 TF
      
      DATA C1 /8.24272D-1/, C2 /2.54627D-3/, C3 /9.848539D-8/
      
      IF (T .GE. 273.15D0) T = T-273.15D0
      TF = T*1.8D0+32.D0
      SVLAMBDA = (C1+C2*TF+C3*TF**2.D0)*1.0D-2
      SVLAMBDA = SVLAMBDA*1.7308D0
      
      END FUNCTION
      

C     ==================================================================       
      REAL*8 FUNCTION SVCP(T)
C     Correlate the specific heat (heat capacity) of saturated stream
C     from "Fundamentals of Classical Thermodynamics - SI Version"
C     by Van Wylen and Sonntag, Table A9, pg. 683
C     Valid for T between 300-3500 K, max error 0.43%
C     T - (input) given temperature [K or C]
C     SVCP - (output) specific heat [J/kg-K]
      REAL*8 T
      REAL*8 C1, C2, C3, C4, E1, E2
      REAL*8 TR
      
      DATA C1 /143.05D0/, C2 /-183.54D0/, C3 /82.751D0/, C4 /-3.6989D0/,
     &     E1 /0.25D0/, E2 /0.5/
     
      IF (T .LT. 273.15D0) T = T+273.15D0
      IF (T .LT. 300.D0) WRITE(*, *) 'Warning: correlation of specific
     &                                heat is out of range'
      TR = T/100.D0
      SVCP = (C1+C2*TR**E1+C3*TR**E2+C4*TR)/18.015D-3
      END FUNCTION            


C     ==================================================================       
      REAL*8 FUNCTION SVP2(T, X)
C     Correlate the vapor pressure [Pa] of aqueous sodium chloride, 
C     by Y. Shibue, Fluid Phase Equilibria, 213(2003) 39-51
C     T - given temperature [K or C]
C     X - given mole fraction of NaCl
      REAL*8 T, X
      REAL*8 TC, PC, TR, GT, HX
      REAL*8 HALF, ONE, TWO, THREE, FOUR, SEVEN
      REAL*8 G1, G2, G3, G4, G5, G6
      REAL*8 A1, A2, A3, B1, B2, U
      
      DATA HALF /.5D0/, ONE /1.D0/, TWO /2.D0/, THREE /3.D0/, 
     &     FOUR /4.D0/, SEVEN /7.D0/
C     Critical properties of water: TC [K] and PC [MPa]
      DATA TC /647.096D0/, PC /22.064D0/
      DATA G1 /-7.85951783D0/, G2 /1.84408259D0/,
     &     G3 /-11.7866497D0/, G4 /22.6807411D0/,
     &     G5 /-15.9618719D0/, G6 /1.80122502D0/
C     Empirial parameters
      DATA A1 /1.28746D-1/, A2 /-7.31097D-1/, A3 /-3.15058D2/
      DATA B1 /3.92767D2/, B2 /-2.46440D3/, U /0.024D0/
      
      IF (T .LT. 273.15D0) T = T-273.15D0
      
      TR = T/TC
      GT = ONE/TR*(G1*(ONE-TR)+G2*(ONE-TR)**(ONE+HALF))+
     &     ONE/TR*(G3*(ONE-TR)**THREE+G4*(ONE-TR)**(THREE+HALF))+
     &     ONE/TR*(G5*(ONE-TR)**FOUR+G6*(ONE-TR)**(SEVEN+HALF))
     
      IF (X .LE. 0.024) THEN
        HX = A2*X/(X+A1**TWO)+A3*X**TWO
      ELSE IF (X .LE. 0.117) THEN
        HX = (A1**TWO*A2/(U+A1**TWO)**TWO+TWO*A3*U)*(X-U)+
     &       B1*(X-U)**TWO+B2*(X-U)*(X**TWO-U**TWO)+
     &       A2*U/(U+A1**TWO)+A3*U**TWO
      ELSE
        WRITE(*, *) 'ERROR: NACL MOLE FRACTION OUT OF RANGE'
        STOP
      END IF
      
      SVP2 = DEXP(DLOG(PC)+GT+HX)*1.D6
      
      END FUNCTION